library(raster)
library(horizon)
library(rgdal)
pro=CRS("+init=epsg:28992")
WGS84<-CRS("+init=epsg:4326")
The altitude data has a resolution of 90 meters.
#r <- getData('alt', country='Austria') #used for downloading the data
r <- raster("data/AUT_msk_alt.grd")
plot(r,main="Altitude map of Austria")
print(r)
## class : RasterLayer
## dimensions : 336, 948, 318528 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 9.4, 17.3, 46.3, 49.1 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : /nobackup/users/dirksen/SVF_highres/SVF/data/AUT_msk_alt.grd
## names : AUT_msk_alt
## values : 108, 3533 (min, max)
#r.NED <- getData('alt', country='NL')# used for downloading the data
r.NED <- raster("data/NLD_msk_alt.grd")
plot(r.NED,main="Altitude map of Netherlands")
print(r.NED)
## class : RasterLayer
## dimensions : 348, 480, 167040 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 3.3, 7.3, 50.7, 53.6 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : /nobackup/users/dirksen/SVF_highres/SVF/data/NLD_msk_alt.grd
## names : NLD_msk_alt
## values : -8, 286 (min, max)
s <- svf(r, nAngles=16, maxDist=500, ll=TRUE)
## 0%...10%...15%...20%...25%...30%...35%...40%...45%...50%...55%...60%...65%...70%...75%...80%...85%...90%...95%...
s.NED <- svf(r.NED, nAngles=16, maxDist=500, ll=TRUE)
## 0%...10%...15%...20%...25%...30%...35%...40%...45%...50%...55%...60%...65%...70%...75%...80%...85%...90%...95%...
print(s)
## class : RasterLayer
## dimensions : 336, 948, 318528 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 9.4, 17.3, 46.3, 49.1 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : /tmp/Rtmp0NECCC/raster/r_tmp_2016-11-24_151201_15395_28883.grd
## names : layer
## values : 0.7161961, 1 (min, max)
print(s.NED)
## class : RasterLayer
## dimensions : 348, 480, 167040 (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : 3.3, 7.3, 50.7, 53.6 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84
## data source : /tmp/Rtmp0NECCC/raster/r_tmp_2016-11-24_151207_15395_49376.grd
## names : layer
## values : 0.9967596, 1 (min, max)
plot(s,main="SVF for Austria")
plot(s.NED,main="SVF for Netherlands")
The raw data from the AHN is available with this link. As an example I included the map number r50fn2.
Besides this product from PDOK and nationaal Georegister also a 3D map of the Netherlands is available
xml.link.raw<-"http://geodata.nationaalgeoregister.nl/ahn2/atom/ahn2_05m_ruw.xml"
r.NED<-raster("AHN/r09dn1.tif") #runnning the 0.5m tif takes approximately 2 hours
plot(r.NED,main="Altitude map of Netherlands")
print(r.NED)
## class : RasterLayer
## dimensions : 12500, 10000, 1.25e+08 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : 110000, 115000, 556250, 562500 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs
## data source : /nobackup/users/dirksen/SVF_highres/SVF/AHN/r09dn1.tif
## names : r09dn1
The runtime for a 0.5m resolution grid is long. Therefore I regrid the raster to a 5m grid. Calculating the SVF will now take approximately 1 minute.
r.NED<-aggregate(r.NED,fact=10) #regridding to a 5m grid, this takes just a minute to run
res(r.NED)
## [1] 5 5
plot(r.NED)
system.time(
s.NED <- svf(r.NED, nAngles=16, maxDist=100, ll=F)
)
## 0%...10%...15%...20%...25%...30%...35%...40%...45%...50%...55%...60%...65%...70%...75%...80%...85%...90%...95%...
## user system elapsed
## 125.513 0.109 126.084
plot(s.NED,main="SVF for Netherlands")
see for the ahn_units.
ahn.units<-readOGR(dsn="AHN_kaartbladen",layer="ahn_units")
## OGR data source with driver: ESRI Shapefile
## Source: "AHN_kaartbladen", layer: "ahn_units"
## with 1441 features
## It has 3 fields
proj4string(ahn.units)<-pro
plot(ahn.units)
coords.knmi<-readRDS("data/coordsKNMI.rda")
coords.knmi<-na.omit(coords.knmi)
coords.knmi<-data.frame(coords.knmi)
coordinates(coords.knmi)<-~DS_LON+DS_LAT
proj4string(coords.knmi)<-WGS84
coords.knmi<-spTransform(coords.knmi,pro)
ov<-over(coords.knmi,ahn.units)
combined<-cbind(ov,as.data.frame(coords.knmi))
print(combined[which(combined$DS_CODE=="260_H"),])
## LO_X LO_Y UNIT DS_CODE DS_LON DS_LAT
## 10 140000 456250 32cn1 260_H 140783.3 456758.2
r.Bilt<-raster("AHN/r32cn1.tif") #runnning the 0.5m tif takes approximately 2 hours
plot(r.Bilt,main="Altitude map of De Bilt")
points(coords.knmi,col="red",lwd=3)
print(r.Bilt)
## class : RasterLayer
## dimensions : 12500, 10000, 1.25e+08 (nrow, ncol, ncell)
## resolution : 0.5, 0.5 (x, y)
## extent : 140000, 145000, 456250, 462500 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725 +units=m +no_defs
## data source : /nobackup/users/dirksen/SVF_highres/SVF/AHN/r32cn1.tif
## names : r32cn1
r.Bilt<-aggregate(r.Bilt,fact=10) #regridding to a 5m grid, this takes just a minute to run
res(r.Bilt)
## [1] 5 5
plot(r.Bilt)
points(coords.knmi,col="red",lwd=3)
system.time(
s.Bilt <- svf(r.Bilt, nAngles=16, maxDist=100, ll=F)
)
## 0%...10%...15%...20%...25%...30%...35%...40%...45%...50%...55%...60%...65%...70%...75%...80%...85%...90%...95%...
## user system elapsed
## 126.430 0.147 127.049
plot(s.Bilt,main="SVF for De Bilt")
points(coords.knmi,col="red",lwd=1)